/****************************************************************/
/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
/*                                                              */
/*          All contents are licensed under LGPL V2.1           */
/*             See LICENSE for full restrictions                */
/****************************************************************/

#include "FreeAnisoHeatConductionMaterial.h"
#include "Function.h"
#include "MooseMesh.h"
#include "Assembly.h"

#include "libmesh/quadrature.h"
registerMooseObject("RoshanApp", FreeAnisoHeatConductionMaterial);
template<>
InputParameters validParams<FreeAnisoHeatConductionMaterial>()
{
	InputParameters params = validParams<Material>();

  params.addCoupledVar("temperature", "Coupled Temperature");

  params.addParam<Real>("thermal_conductivity_1", "The thermal conductivity in the 1st main direction");
  params.addParam<Real>("thermal_conductivity_2", "The thermal conductivity in the 2nd direction");
  params.addParam<Real>("thermal_conductivity_3", "The thermal conductivity in the 3rd direction");
  params.addParam<Point>("direction_1", "The vector of the 1st main direction");
  params.addParam<Point>("direction_2", "The vector of the 2nd main direction");
  params.addParam<Point>("direction_3", "The vector of the 3rd main direction");
  params.addParam<Real>("specific_heat", "The specific heat value");
  params.addParam<Real>("density", "The _density value");


  return params;
}



FreeAnisoHeatConductionMaterial::FreeAnisoHeatConductionMaterial(const InputParameters & parameters) :
  Material(parameters),

  _has_temperature(isCoupled("temperature")),
  _temperature(_has_temperature ? coupledValue("temperature") : _zero),
  _my_thermal_conductivity_1(getParam<Real>("thermal_conductivity_1")),
  _my_thermal_conductivity_2(getParam<Real>("thermal_conductivity_2")),
  _my_thermal_conductivity_3(getParam<Real>("thermal_conductivity_3")),
  _my_direction_1(getParam<Point>("direction_1")),
  _my_direction_2(getParam<Point>("direction_2")),
  _my_direction_3(getParam<Point>("direction_3")),

  _my_specific_heat(isParamValid("specific_heat") ? getParam<Real>("specific_heat") : 0),
  _my_density(isParamValid("density") ? getParam<Real>("density") : 1),



//  _translate_conductivity(declareProperty<RankTwoTensor>("translate_conductivity")),

  _direction_1(declareProperty<Point>("direction_1")),
  _direction_2(declareProperty<Point>("direction_2")),
  _direction_3(declareProperty<Point>("direction_3")),
  _thermal_conductivity_1(declareProperty<Real>("thermal_conductivity_1")),
  _thermal_conductivity_2(declareProperty<Real>("thermal_conductivity_2")),
  _thermal_conductivity_3(declareProperty<Real>("thermal_conductivity_3")),

//  _translate_conductivity(declareProperty<RankTwoTensor>("translate_conductivity")),
//  _conductivity_matrix(declareProperty<RankTwoTensor>("conductivity_matrix")),
//  _translate_conductivity_Trans(declareProperty<RankTwoTensor>("translate_conductivity_Trans")),


  _specific_heat(declareProperty<Real>("specific_heat")),
  _density(declareProperty<Real>("density")),
_conductivity_matrix_xyz(declareProperty<RankTwoTensor>("thermal_conductivity_xyz"))
{

}

Real
FreeAnisoHeatConductionMaterial::compute_cos(Point p1,Point p2)
{
	Real cos_value;
	cos_value = p1*p2/p1.norm()/p2.norm();
	return (cos_value);
}

//RankTwoTensor
//FreeAnisoHeatConductionMaterial::translate_conductivity(Real S11, Real S21, Real S31, Real S12, Real S22, Real S32, Real S13, Real S23, Real S33)
//{
//	RankTwoTensor _translate_conductivity;
//	 _translate_conductivity.RankTwoTensor(S11, S21, S31, S12, S22, S32, S13, S23, S33);
//	return _translate_conductivity;
//}
//
//RankTwoTensor
//FreeAnisoHeatConductionMaterial::conductivity_matrix(Real S11, Real S21, Real S31, Real S12, Real S22, Real S32, Real S13, Real S23, Real S33)
//{
//	RankTwoTensor conductivity_matrix;
//	conductivity_matrix.RankTwoTensor(S11, S21, S31, S12, S22, S32, S13, S23, S33);
//	return conductivity_matrix;
//}




void
FreeAnisoHeatConductionMaterial::computeProperties()
{
  for (unsigned int qp(0); qp < _qrule->n_points(); ++qp)
 {
   _thermal_conductivity_1[qp] = _my_thermal_conductivity_1;
   _thermal_conductivity_2[qp] = _my_thermal_conductivity_2;
   _thermal_conductivity_3[qp] = _my_thermal_conductivity_3;


   _direction_1[qp] = _my_direction_1;
   _direction_2[qp] = _my_direction_2;
   _direction_3[qp] = _my_direction_3;
   Point nx(1,0,0);
   Point ny(0,1,0);
   Point nz(0,0,1);

   Real l1 = compute_cos(_direction_1[qp],nx);//l1 m1 n1 见文献《正交各向异性材料三维热传导问题的有限元列式》，闫相桥等
   Real l2 = compute_cos(_direction_2[qp],nx);
   Real l3 = compute_cos(_direction_3[qp],nx);
   Real m1 = compute_cos(_direction_1[qp],ny);
   Real m2 = compute_cos(_direction_2[qp],ny);
   Real m3 = compute_cos(_direction_3[qp],ny);
   Real n1 = compute_cos(_direction_1[qp],nz);
   Real n2 = compute_cos(_direction_2[qp],nz);
   Real n3 = compute_cos(_direction_3[qp],nz);
//   std::cout<<l1<<" "<<l2<<" "<<l3<<" "<<m1<<" "<<m2<<" "<<m3<<" "<<n1<<" "<<n2<<" "<<n3<<std::endl;
//   RankTwoTensor A(_thermal_conductivity_1[qp], _thermal_conductivity_1[qp], _thermal_conductivity_1[qp]);
   RankTwoTensor translate_conductivity(l1,m1,n1,l2,m2,n2,l3,m3,n3);
   RankTwoTensor conductivity_matrix(_thermal_conductivity_1[qp],0,0,0,_thermal_conductivity_2[qp],0,0,0,_thermal_conductivity_3[qp]);
   RankTwoTensor translate_conductivity_Ttrans(l1,l2,l3,m1,m2,m3,n1,n2,n3);
   RankTwoTensor conductivity_matrix_xyz = translate_conductivity*conductivity_matrix*translate_conductivity_Ttrans;
   _conductivity_matrix_xyz[qp] = conductivity_matrix_xyz;
//   std::cout<<_conductivity_matrix_xyz[qp](0,0)<<","<<_conductivity_matrix_xyz[qp](0,1)<<","<<_conductivity_matrix_xyz[qp](0,2)<<std::endl;
//   std::cout<<_conductivity_matrix_xyz[qp](1,0)<<","<<_conductivity_matrix_xyz[qp](1,1)<<","<<_conductivity_matrix_xyz[qp](1,2)<<std::endl;
//   std::cout<<_conductivity_matrix_xyz[qp](2,0)<<","<<_conductivity_matrix_xyz[qp](2,1)<<","<<_conductivity_matrix_xyz[qp](2,2)<<std::endl;
   _specific_heat[qp] = _my_specific_heat;
   _density[qp] = _my_density;
 }
}

